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ABSTRACT 

We present observations and dynamical models of the stellar nuclear clusters (NCs) at the 
centres of NGC 4244 and M33. We then compare these to an extensive set of simulations test- 
ing the importance of purely stellar dynamical mergers on the formation and growth of NCs. 
Mergers of star clusters are able to produce a wide variety of observed properties, including 
densities, structural scaling relations, shapes (including the presence of young discs) and even 
rapid rotation. Nonetheless, difficulties remain, most notably that the second order kinematic 
moment Vmis = Vv^ + ct^ of the models is too centrally peaked to match observations. 
This can be remedied by the merger of star clusters onto a pre-existing nuclear disc, but the 
line-of-sight velocity V is still more slowly rising than in NGC 4244. Our results therefore 
suggest that purely stellar dynamical mergers cannot form NCs, and that gas dissipation is 
a necessary ingredient for at least ^ 50% of a NCs mass. The negative vertical anisotropy 
found in NGC 4244 however requires at least 10% of the mass to be accreted as stars, since 
gas dissipation and in situ star formation leads to positive vertical anisotropy. 
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1 INTRODUCTION 

Studies of the centres of galaxies across the Hubble sequence have 
shown that they frequently host central massive objects (CMOs) 
such as supermassive black holes (SMBHs) and massive stellar nu- 
clear clusters (NCs). NCs are very common: the ACS Virgo Clus- 
ter Survey found a NC in 66 — 88% of low and intermediate lu- 
minosity early-type galaxies (Grant et al. 2005; Cote et al. 2006). 
A similar fraction, ~ 75%, of late-type spiral galaxies contain a 
NC (Boker et al. 2002). However NCs do not appear to populate 
the most massive early-type galaxies brighter than Mb ~ —20.5 
(Grant et al. 2005; Cote et al. 2006). Both Rossa et al. (2006) and 
Cote et al. (2006) found that the luminosity of the NC correlates 
with that of the host galaxy in elliptical and early-type spiral galax- 
ies. Remarkably, both NCs and SMBHs follow the same correlation 
between their mass, Mcmo , and that of their host galaxy, with NCs 
on the lower end and SMBHs on the upper end of the relation (Fer- 



E-mail: mhartmannSuclan . ac . uk 
t E-mail: vpdebattistaSuclan . ac . uk, RCUK Fellow 
I E-mail: asethScf a . harvard . edu, CfA Fellow 
§ E-mail: cappellari@astro.ox.ac.uk 
^ E-mail: trq@astro . Washington . edu 



rarese et al. 2006b; Wehner & Harris 2006). Likewise, both NC 
and SMBH masses follow an Mcmo — ""e relation with the veloc- 
ity dispersion of their host spheroid (Ferrarese et al. 2006b; Wehner 
& Harris 2006), with the relation for NCs parallel to, but ~ lOx 
more massive at a given CTe, that for SMBHs. These observational 
facts suggest that SMBHs and NCs share a similar growth history 
(McLaughlin et al. 2006). This motivates us to understand SMBH 
growth by using the formation of NCs as their visible proxies. An 
important advantage of NCs is that there are additional observables 
which can be measured for them, in addition to their mass. In par- 
ticular one can characterise their kinematics and stellar population. 
This provides further constraints to their formation, which we try 
to exploit in this paper. 

Two main scenarios have been proposed to explain the for- 
mation of NCs. One scenario relies on NCs forming in situ 
out of gas falling into the centre. A number of mechanisms 
have been proposed for driving gas to galactic centres, including 
the magneto-rotational instability (Milosavljevic 2004), gas cloud 
mergers (Bekki 2007) or the action of instabilities (Shlosman & 
Begelman 1989; Maciejewski et al. 2002; Schinnerer et al. 2003, 
2008). Alternatively, NCs may form from the merging of star clus- 
ters (SCs) after sinking to the centre under the action of dynamical 
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friction (Tremaine et al. 1975; Miocchi et al. 2006). Indeed Lotz 
et al. (2001) found a depletion of bright globular clusters within the 
inner region of dEs and similar colours for the globular cluster and 
NC population which, however, are bluer than the underlying stel- 
lar population of the host. Analysis of NC colours in dE galaxies by 
Lotz et al. (2004) and Cote et al. (2006) suggest that many, but not 
all dE nuclei could be explained by a cluster merger scenario. Self- 
consistent simulations have found that mergers can occur and the 
resulting NC masses and sizes are consistent with those observed 
(Bekki et al. 2004; Capuzzo-Dolcetta & Miocchi 2008a,b). 

NC formation in late-type galaxies is an ongoing process, with 
star formation histories that are extended, possibly constant (Rossa 
et al. 2006; Walcheretal. 2006; Scth ct al. 2010), and stars younger 
than 100 Myrs present. Mergers and accretion of old SCs is thus not 
a viable formation scenario, regardless of whether enough such SCs 
are available. 

However, mergers and accretion of young SCs formed near 
the centres of galaxies, such as those observed in the Milky Way 
(Figer et al. 1999, 2002), NGC 2139 (Andersen et al. 2008), and 
NGC 253 (Komei & McCrady 2009), could still provide a viable 
formation mechanism. Agarwal & Milosavljevic (201 1) used an- 
alytic modelling to show that infalling SCs from an empirical SC 
population produce NCs of the right mass in isolated spheroidal and 
late-type galaxies. Such SCs must form quite close to the galaxy 
centres (<1 kpc), otherwise the timescales for infall due to dynam- 
ical friction are prohibitively long (Milosavljevic 2004). In a sam- 
ple of observed bulgeless spiral galaxies Neumayer et al. (2011) 
find that the dynamical friction timescales are < 2 Gyrs for SCs 
with masses > 2 x 10^ M© within 500 pc. The supply of young 
SCs to the inner regions may be enhanced because tidal forces are 
compressive within 10% of the effective radius when the density 
profile has a Sersic index n < 2 (Emsellem & van de Ven 2008). 
However, a galaxy with a constant density dark matter core would 
inhibit infall altogether (Read et al. 2006; Goerdt et al. 2008, 2010). 

The morphology and kinematics of NCs in disc galaxies pro- 
vide additional constraints on their formation. From observations 
of edge-on galaxies, Seth et al. (2006) find that their NCs are typ- 
ically elongated along the plane of the galaxy disc. These NCs are 
often compound structures, with a younger, bluer thin disc (NCD) 
embedded within an older spheroidal component (NCS). In the case 
of the NC in NGC 4244, spectra reveal young (< 100 Myr) stars in 
the NCD amounting to ~ 5% of the total NC mass. Furthermore, 
integral-field spectroscopy of this NC shows rotation in the same 
sense as the galaxy (Seth et al. 2008b). As we show below, it is not 
yet possible to determine whether this rotation is restricted to the 
NCD. 

Accretion without gas dissipation need not result in slowly 
rotating systems: Read et al. (2008) showed that satellite galaxies 
that accrete onto disc galaxies can get dragged down on to the plane 
of the disc and settle directly into a rapidly rotating thick disc with 
(V/ct), as large as 4. Likewise, Eliche-Moral et al. (2006) find that 
low density satellites can accrete onto bulges and end up rapidly 
rotating. Can NCDs form through an analogous process? Since SCs 
falling into the nucleus are likely to have formed in the mid-plane 
of the galaxy, they will generally share in its rotation. 

In this Paper we explore in detail whether the SC 
merger/accretion scenario is able to produce the density and kine- 
matic properties observed in NCs. In Section 2 we present the ob- 
servations of the NCs in NGC 4244 and M33. In Section 3 we 
model integral field data of NGC 4244 and M33 using the Jeans 
Anisotropic MGE method of Cappellari (2008). In Section 4 we 
describe our simulations and methods. We present the results of 



NCs formed by mergers/accretions of SCs in Section 5. In Section 
6 we introduce an initial NCD and follow its evolution as it accretes 
SCs. In Section 7 we apply the Jeans Anisotropic MGE modelUng 
also to the simulations. We discuss the implications of our simula- 
tions and show that a purely collisionless merger origin of NCs is 
not consistent with the observations in Section 8. 



2 OBSERVATIONAL CONSTRAINTS 

We compare our simulations with observations of the NCs in two 
nearby galaxies, M33 and NGC 4244. These two galaxies pro- 
vide us with well-resolved nuclei for which integral field kinematic 
data are available. While M33 is inclined by i — 49° (Corbelli & 
Schneider 1997), NGC 4244 is an edge-on galaxy. These two ob- 
jects thus give us two different perspectives for studying the mor- 
phology and kinematics of NCs in disc galaxies. In this paper we 
assume a distance of 4.4 Mpc for NGC 4244 (Seth et al. 2005) and 
0.8 Mpc for M33 (Lauer et al. 1998). 

2.1 Spectroscopy of NGC 4244 and M33 

Seth et al. (2006) presented Hubble Space Telescope F814W pho- 
tometry of the NC in NGC 4244, which is resolved into a nuclear 
cluster spheroid (NCS) and a bluer nuclear cluster disc (NCD). The 
half mass radius, -Rcff , of the NC obtained by fitting a King profile 
to the ACS/F814W images (Seth et al. 2006) is about 5 pc. 

Here we use K-hmd spectra of NGC 4244 (Seth et al. 2008b) 
and M33 (Seth et al, in prep) obtained at Gemini North with the 
Near-Infrared Integral Field spectrograph (NIFS), an image slic- 
ing field unit spectrograph. The PSF core in both observations is 
~ 0.1"FWHM, with final data cubes sampled with 0'.'05 pixels. In 
NGC 4244, this corresponds to 1.06 pc pixel"^, while in M33 the 
pixel scale is 0.19 pc pixel" ^. The imperfect adaptive optics cor- 
rections result in nearly half of the light being scattered PSF halo 
0'.'7 FWHM). This leads to problems separating the rotation in 
the NCS and NCD in NGC 4244 (see Section 3.2). In both clusters, 
the rotation amplitude increases with radius out to the effective ra- 
dius, and decreases at larger radii. However, whether or not a NCD 
is present in the NC of M33 cannot be determined because of the 
galaxy's inclination. In Figure 1 we show the line-of-sight velocity 
profiles of M33 and NGC 4244. 

2.2 Isophotal shape of the NC in NGC 4244 

The shape of an isophote can be quantified by the Fourier coeffi- 
cients of the expansion 

I{(j)) = Io + Y^ An sin(n0) + B„ cos(n0) (1) 

n=l 

where 7o is the mean intensity along the ellipse, </> is the az- 
imuthal angle and An, B„ (n = 1, 2, ...) are harmonic ampli- 
tudes. The ellipse which best fits the isophote has the coefficients 
An, Bn (n — 1,2) equal to zero. The deviations from the best- 
fit isophote are then given by the higher order coefficients An, 
B„ (n = 3, 4, ...). The leading residual term generally will be the 
n = 4 term, which determines whether the isophote is discy 
(B4 > 0) or boxy {B4 < 0). We use the task ELLIPSE in IRAF to 
measure B4. The parameter B4 has been shown to correlate with 
kinematic properties of the host galaxy in observations (Bender 
1988; Kormendy & Djorgovski 1989; Kormendy & Bender 1996) 
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Figure 1. Profiles from the NIFS data (Seth et al. 2()()8b) of line-of-sight 
velocities along the major axis in NGC 4244 (top) and M33. The kinematics 
were extracted from the integral-field data within a slit of |z| < 0.4 pc. This 
corresponds to 0.08 Rgg for NGC 4244 and 0.2 R^g for M33. 
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Figure 2. B4 for the NC in NGC 4244 derived from the integrated K-hmd 
flux data (Seth et al. 2008b). 



and in simulations of galaxy mergers (Naab et al. 1999; Naab & 
Burkert 2003; Naab et al. 2006). This result was strongly confirmed 
with the large ATLAS^° sample of early-type galaxies (Cappellari 
et al. 201 1) which demonstrated that a significant disciness is an 
unambiguous indication of a fast rotator galaxy (Emsellem et al. 
201 1). Figure 2 shows B4, for the NC in NGC 4244, which is clearly 
discy. 



2.3 Axisymmetry of the M33 nucleus 

The face-on axial S3fmmetry is an important constraint on the for- 
mation of any stellar system. In general the merger of stellar sys- 
tems leads to triaxial systems, e.g. elliptical galaxies are found to 
be triaxial (Wagner et al. 1988; Naab et al. 1999; Naab & Burkert 
2003; Naab et al. 2006). Very Uttle data are available for the face- 
on axial symmetry of NCs. For the NC in M33 we find evidence 
for axisymmtry. The inner 0.5" (= 1.8 pc) of the NC in M33 has 
an apparent b/a = 0.84, with an average PA = 19.3° (Lauer et al. 
1998), which is close to the PA of the inner disc PA = 23° ± 1° 
within R < 4.0 kpc (Zaritsky et al. 1989). Thus we can assume that 
the NC is in the same plane as the inner disk and therefore it is pos- 
sible to determine the NCs apparent axis ratio b/a. The galaxy is 
inclined by i = 49° (Corbelli & Schneider 1997), so that an oblate 
spheroid of apparent axis ratio b/a at an inclination i would have 
an intrinsic vertical-to-horizontal axis-ratio go given by 



go = 



{b/aY — cos^ 



(Hubble 1926). The go ~ 0.7 that Eqn. 2 implies for M33's NC 
is fully consistent with the range of values of go found in edge- 
on galaxies by Seth et al. (2006) (NGC 4244 having go ~ 0.5). 
Moreover we measured the kinematical PAkin = 36 ± 18 (Icr 
error) from the NIFS integral-field kinematics, with the routine 
FIT_KINEMATIC_PA described in Appendix C of Krajnovic et al. 
(2006). Thus the kinematical misalignment is consistent with zero, 
consistent with the NC of M33 being axisymmetric. 



3 MODELLING THE OBSERVATIONS 

We want to compare the kinematics of the observed NCs to those of 
the simulations. Given that integral-field kinematics are available 
for the observed NCs, we use two methods which make full use 
of the two-dimensional data. One method is based on the Jeans 
equations and the other is based on the (V/u, e) diagram. In the 
following we show that consistent results are obtained with both 
approaches. 



3.1 The Jeans Anisotropic MGE dynamical models 

The integral-field stellar kinematics for the NCs in NGC 404 (Seth 
et al. 2010), NGC 4244 (Seth et al. 2006) and M33 (Seth et al, in 
prep), the few NCs for which this type of data are available, suggest 
that NCs are likely not far from oblate axisymmetric systems (as we 
also argued in Section 2.3). A classic reference model to quantify 
the rotation of axisymmetric galaxies is the isotropic rotator (Bin- 
ney 1978), to which real galaxies have often been compared (e.g. 
Satoh 1980; Binney et al. 1990; van der Marel 1991 ; van der Marel 
& van Dokkum 2007; van der Wei & van der Marel 2008). 

Recent dynamical modelling studies have found that real ax- 
isjonmetric galaxies are generally best matched by models which 
have a nearly oblate velocity ellipsoid, flattened along the direc- 



tion of the symmetry axis ^ on 



(Cappellari et al. 



(2) 



2007; Cappellari 2008; Thomas et al. 2009). A useful general- 
isation of the isotropic rotator to quantify the rotation of these 
anisotropic systems is then a rotator with oblate velocity ellipsoid 
(az ^ (JR = CT^), which provides a good approximation for the ob- 
served integral-field kinematics of real galaxies (Cappellari 2008; 
Scott et al. 2009). 

To perform the measurement of the degree of rotation we used 
the Jeans Anisotropic MGE (JAM) software^ which implements an 
efficient solution of the Jeans equations with oblate velocity ellip- 
soid (Cappellari 2008). Under that assumption the model gives a 
unique prediction for the observed first two velocity moments V 
and Vrms = VV^ -t- u^, where V is the observed mean stellar ve- 
locity and cr is the mean stellar velocity dispersion. 

The luminous matter likely doininates in the high-density nu- 
clei of the studied galaxies. The same is true by construction for the 
simulations. For this reason the dynamical models assume that light 
traces mass. To parametrise the surface brightness distribution of 
either the galaxies or the iV-body simulations we adopted a Multi- 
Gaussian Expansion (MGE; Emsellein et al. 1994), which we fit 
to the images with the method and software^ of Cappellari (2002). 
For a given inclination the models then have one free nonlinear 
parameter, the anisotropy fiz — 1 — ai/a\, and two linear scal- 
ing factors: (i) the dynamical M/ L and (ii) the amount of rotation 
K = Lobs/iobi, which is the ratio between the observed projected 
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angular momentum Lobs and the one for a model with oblate veloc- 
ity ellipsoid Lobi (see Cappellari 2008, for details). To find the best 
fitting model parameters we constructed a grid in the non-Unear pa- 
rameter /3z and for each value we linearly scaled the M/L to match 
the Vrms data in a 'x^ sense. At the best-fitting {j3z, M/L) we then 
computed the model velocity V, further assuming or = o<t>, and 
measure the amount of rotation k from the observed velocity. 

3.2 JAM models of observed NCs 

We applied the procedure to model the Gemini NIFS stellar kine- 
matics of the NC in NGC 4244 (Seth et al. 2008b). The MGE fit 
to the nuclear part of the ACS/F814W image is shown in Figure 3. 
We used the ACS/F814W image for the MGE fit, because of the 
better known point spread fimction compared to the /f-band NIFS 
data. The NC was assumed to be a dynamically distinct component, 
in equilibritmi in the combined potential due to the galaxy and the 
NC itself We used all MGE Gaussians, of both the NC and the 
main galaxy disc, in the calculation of the gravitational potential, 
but only the nuclear Gaussians were used to parametrise the sur- 
face brightness of the NC. Thus our JAM models of NCs are not 
self-consistent. Although this assumption is physically motivated 
and can sometimes produce significant differences in the kinemat- 
ics of the model, in this case the results are indistinguishable from a 
self-consistent model. The best-fitting JAM model for NGC 4244, 
which has an edge-on inclination (i = 90°), is shown in Figure 4. 
It has a best-fitting anisotropy parameter Pz = — 0.2±0.1 and a 
rotation parameter n — 0.99 ± 0.05. This implies that the NC ro- 
tates almost perfectly as the reference rotator with oblate velocity 
ellipsoid (for which k = 1). The best-fitting JAM model has no 
central intermediate mass black hole (IMBH). However from our 
NIFS data we can place an upper limit of M. < 1 x 10 ' Mq on 
the mass of a possible IMBH inside the NC. This is < 1% of the 
mass of the cluster, Mnc = (1-1 ± 0.2) x 10^ M©. For larger 
M, the model shows a clear central peak in Vrms which is strongly 
excluded by the data. 

The NC of NGC 4244 is composed of a NCD and a NCS (Seth 
et al. 2006), which have distinct stellar populations (Seth et al. 
2008b). While Seth et al. (2008b) argued that both the disc and 
spheroid are rotating due to the change of line strengths and optical 
colour, the models suggest that the rotation signal may come only 
from the disc. For this we fitted two MGE models separately for 
the spheroidal and disc component defined by Seth et al. (2006). 
We computed the best-fitting JAM model for the Vrms like before, 
but then we fit V by setting the rotation to zero (k — 0) for the 
MGE Gaussians describing the NCS. This model with an unro- 
tating NCS reproduces the NIFS data as well as the model with 
constant anisotropy and rotation for both the NCS and the NCD. 
This is because the flat disc component of the NC dominates the 
light and the rotation of the model in the region where NIFS data 
are available. More spatially resolved/extended kinematics would 
be required to measure the rotation of the spheroidal component of 
theNC. 

We constructed a similar MGE model for the NC of M33 using 
the WFPC2/F814W photometry. For the JAM model we adopted 
the inclination of the main galaxy disc (i « 49°; Corbelli & 
Schneider 1997). In this case the observed distribution of stellar a 
from the NIFS data is quite irregular and presents significant asym- 
metries, which cannot be reproduced by an axisymmetric model. 
The irregularity in the velocity field is possibly due to granular- 
ity in the velocity field, with individual bright AGB stars having 
significant influence on the velocity and dispersion measurements. 
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Figure 3. The contours of the surface brightness of the ACS/F814W image 
of the NC in NGC 4244, in steps of 0.5 mag, are overlaid on the PSF- 
convolved MGE model. Both the NCD and the NCS are well described by 
the model. The scale is ~ 21.1 pc/arcsec with a half-mass radius ~ 0.27" 
(Seth et al. 2008b). The dashed box shows the region within which kine- 
matic data are observed. 



The velocity field of M33 and other nearby NCs will be discussed 
in more detail in an upcoming paper (Seth et al. in prep). As the 
anisotropy cannot be reliably inferred from the data, we assumed 
isotropy and derived the dynamical AI/L from a fit to the observed 
Vrms. The rotation parameter derived from the observed V is then 
K = 1.02 ± 0.10, which, as for NGC 4244, is consistent with the 
rotation of the reference rotator with oblate velocity elUpsoid. The 
total mass of the NC is Mnc = (1-4 ± 0.2) x 10*^ Mq. 

3.3 Rotation from the (V/a, e) diagram 

An alternative classic way of quantifying rotation is given by the 
{V/a, e) diagram (lUingworth 1977; Binney 1978). Traditionally 
the observed V/a quantity was computed from the central veloc- 
ity dispersion and the maximum rotational velocity. Binney (2005) 
updated and improved the formalism to compute the quantity in a 
more robust way when integral-field data are available. Here the 
availability of integral-field kinematics for both the observations 
and the simulations allow us to apply this improved method. We 
use the updated formulae and define 



V 



(3) 



which we applied within 2 Rcu as a luminosity-weighted quan- 
tity, which we estimate from the integral-field kinematics. Here F„ 
is the flux contained inside the n-th Voronoi bin and V„ and a„ 
the corresponding measured mean line-of-sight velocity and veloc- 
ity dispersion. The ellipticity is defined following Cappellari et al. 
(2007) by a similar expression as 



(1 



2 

■ 1 = 



F r2 ' 



(4) 
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Figure 4. JAM model for the stellar kinematics of NGC 4244. The top two panels show the bi-symmetrized NIFS Vrma (left) and V (right). The two 
bottom panels show the corresponding kinematics of the best fitting JAM model. The contours of the surface brightness are overlaid in steps of 0.5 mag. The 
dots indicate the position of the centroids of the Voronoi bins for which the kinematics were extracted. The best fitting model has Pz = —0.2 ±0.1 and 
K = 0.99 ± 0.05. The NC of NGC 4244 has an almost perfect oblate velocity elUpsoid. 



where the [x, y) coordinates are centred on the galaxy nucleus and 
the X axis is aligned with the NC photometric major axis. We esti- 
mate e from the individual pixels, inside a given galaxy isophote, 
within the same region used for the computation of (V/ct)^. The 
main disadvantage of the (V/ct, e) diagram, with respect to the 
JAM models, is that it does not rigorously take into account of mul- 
tiple photometric systems, like a disc and a spheroid. Moreover, 
while the diagram can quantify the anisotropy, it does not provide 
any information on whether it is mostly radial or tangential. Still 
the diagram provides an important independent test of more de- 
tailed models and provides an easy way of comparing simulations 
and observations. 



We used Eqns. 3 and 4 to place the NC of M33 and NGC 4244 
on the (V/(T, e) diagram, using our NIFS data. Given that the dia- 
gram is defined for edge-on orientations, while M33 has an incli- 
nation of i « 49°, we projected the (V/cr, e) values for the NC 
to an edge-on view following Binney & Tremaine (2008). The NC 
in NGC 4244 is seen at a nearly edge-on orientation and is weakly 



anisotropic. The location of M33 is slightly more uncertain, given 
the non edge-on view, but the NC is consistent with isotropy. 



4 NUMERICAL METHODS 

The half mass relaxation time for the NC in NGC 4244, with 
i?eff ~ 5 pc and total mass 1.1 x lO'^ M© is ~ 10 Gyrs. NCs 
in the Virgo Cluster Survey have relaxation times ranging from 
1 — 10 Gyrs (Merritt 2009). Therefore it is reasonable to approx- 
imate NCs as coUisionless systems on timescales of < 1 Gyr, al- 
lowing us to use standard coUisionless codes to simulate their evo- 
lution. 

In this section we describe the A*'-body initial conditions for 
the simulations. Because we are interested in evolution at the inner 
~ 100 pc region of such galaxies, we neglect the dark matter halo. 
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Figure 5. Real and simulated NCs on the (V/a,e) diagram of Biimey 
(2005). The thick green line indicates the location of edge-on isotropic mod- 
els, while the other solid lines are anisotropic models with global anisotropy 
(5 = 1 — 2cr^ / ^cr^ + , spaced at 0. 1 intervals. The magenta line is the 
lower envelope for fast-rotating galaxies defined in Cappellari et al. (2007). 
The dashed lines indicate how the magenta line transforms at lower in- 
clinations in 10° steps, indicated by the dotted lines. The green stars and 
red triangles indicate the location of the simulated NCs, while the magenta 
squares are the observed NCs. The observed location for M33 has been pro- 
jected (red line) assuming an inclination of i = 49°. For comparison we 
also plot the fast-rotator (blue dots) and slow-rotator (red circles) eai'ly-type 
galaxies from Cappellaii et al. (2007). Both the simulated and real NCs are 
fast rotating and relatively close to the isotropic line. 



4.1 Galactic disc model 

NGC 4244 and M33 are late-type, bulgeless galaxies. We consider 
only the main galactic disc and adopt an exponential profile: 



p{R,z) 



-R/Rd 



-.'12.1 



2'KZd 



(5) 



where Md is the disc mass, is the scale-length of the disc and 
Zd the scale-height. We use Ma = 10^° Mq, Rd = 2.0 kpc, 
Zd = 200 pc and truncate the disc at 5 Rd- The disc is repre- 
sented by 4 X 10® particles. If we had used equal mass particles 
for the galactic disc, this would correspond to particle masses of 
2500 Mq. Such high masses would inhibit dynamical friction on 
objects of mass ~ 10* Mq, and lead to excessive heating of any in- 
falling clusters. In order to reduce such effects, we use multi-mass 
disc particles, with masses ranging from 7 Mq within the inner 
20 pc increasing to 1.2 x 10^ Mq in the disc's outskirts. Figure 6 
shows the distributions of masses and softenings of disc particles; 

1/3 

the latter is related to particle mass via oc nip . We use the 
epicyclic approximation to set Toomre-Q = 1.2. This setup is im- 
perfect and needs to be relaxed; moreover simulations using such 
initial conditions can only be run for a few crossing times of the 
main disc, since the radial migration of stars induced by disc in- 
stabilities, including bars and spirals, would introduce massive par- 
ticles to the central regions (Sellwood & Binney 2002; Debattista 
et al. 2006; Roskar et al. 2008). In order to check that the system 
does not rapidly homogenise, we relaxed the system for 10 Myrs, 




10' 10^ 10^ 10* 10^ 10* 10^ 




6„ [pc] 



Figure 6. The uni'elaxed initial conditions of the main galactic disc. Left: 
The solid line shows the number of particles (left axis) and the dashed line 
the radius R (right axis) versus the particle mass. Right: The solid line rep- 
resents the number of particles with a certain softening e and the dashed 
line the softening of the particles versus the radius. 
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Figure 7. The masses of particles in the inner 100 pc of the main disc in 
run M2. Blue: initial setup; black: initial conditions after relaxing the disc 
for 10 Myr; red: after 17.5 Myr; green: after 35 Myr. The top panel shows 
the number of particles of a given mass while the bottom panel shows the 
cumulative mass distribution. 



after which main disc particles in the central 100 pc have masses 
ranging from 7 Mq to 2.8 x 10'^ Mq . In Figure 7 we plot the mass 
distribution for the disc particles within the inner 100 pc at three 
different times for the longest merger simulation (M2, described 
below). The cumulative distribution (bottom panel) shows that the 
mass at the centre is dominated by particles with masses smaller 
than 2500 Mq, with more than half the mass in this region coming 
from particles with masses less than 1000 Mq. 



4.2 Bulge model 

On timescales comparable to the crossing time of the galaxy, dy- 
namical instabilities, such as bars and spirals, move particles to the 
central regions. Multi-mass particle simulations within discs there- 
fore are no longer possible in these cases. We therefore set up a 
bulge, again neglecting the dark matter halo, and evolve the sys- 
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Figure 8. The masses of particles in the inner 30 pc of the bulge in simula- 
tion Al. Black: initial setup; red: initial conditions after the 65 Myr needed 
for the infall of the NCS seed; green: after 0.275 Gyrs. The top panel shows 
the number of particles of a given mass while the bottom panel shows the 
cumulative mass distribution. 



terns in this environment. The bulge model has a Hemquist (1990) 
profile: 



p{r) = 



2-Kr (r + of' 



(6) 



where Mb is the bulge mass and a is the scale radius. We use 
Mi, = 5 X 10^ M0 and a = 1.7 kpc. The bulge is truncated 
by eliminating all particles with enough energy to reach r > 15a, 
therefore the density drops gently to zero at this radius (Sellwood & 
Debattista 2009). The bulge is populated by 3.5 x lO'' particles with 
masses ranging from 40 M© to 3.9 x 10^ M©; particle masses are 
selected by the weighting function w{L) oc 3 + 50001/^, with L the 
specific angular momentum, ensuring a high resolution within the 
inner 160 pc (Sellwood 2008). The softening is related to the parti- 

1 /3 

cle mass via ep oc nip . Unlike the disc, the bulge has no strong in- 
stabilities, therefore the distribution of particles remains unchanged 
even on timescales of a few Gyrs. We need these timescales to 
model multiple accretions of SCs. Figure 8 shows the particles dis- 
tribution for simulation Al, described below, within 32 pc for three 
different times. This simulation shows the largest changes in the 
mass distribution. Although accretion delivers higher mass parti- 
cles into the central region, the distribution of particles does not 
change substantially. 



Table 1. The SCs used in the simulations. M» is the stellar mass of the SC, R^ff 
is the effective (half-mass) radius, c is the concentration (see text for definition) 
and M, is the mass of the black hole if one is present. 



Model 


M, 




c 


M. 




[Mq] 


[pc] 




[M0] 


CI 


2.2 x 10« 


1.72 


0.07 




C2 


2.2 X 10^ 


1.60 


0.07 


4.4 X 10" 


C3 


2.0 X 10'^ 


2.18 


0.12 




C4 


2.0 X 10^ 


1.11 


0.12 




C5 


6.0 X 10^ 


1.11 


0.16 





ro O 
I ^ 
O 
D. 

o 
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Figure 9. Volume density of the SC models CI - C5. 



C3-C4 and 15 M© for C5) and equal softening (e = 0.104 pc for 
CI and C2 and e = 0.13 pc for C3-C5). IMBHs may be present in 
some SCs (Gebhardt et al. 2000, 2005; Gerssen et al. 2002, 2003; 
Noyola et al. 2008, 2010, but see also van der Marel & Anderson 
2010). In model C2 we include an IMBH at the centre of the cluster 
CI by adiabatically growing the mass of a single particle with soft- 
ening Ep = 0.042 pc over 100 Myrs. Table 1 lists the properties of 
the SC models. The concentration c is defined as c = log(_Reft /^c) 
where R^a is the half mass radius (effective radius) and Rc is the 
core radius, where the surface density drops to half of the central. 
Figure 9 plots the volume density profiles of the SCs. The central 
density po ranges from 3 x 10'^ to 1 x 10^ Mq pc~^; the masses 
and Rett (see Tab. 1) are comparable to young massive star clusters 
in the Milky Way (Figer et al. 1999, 2002), in the LMC (Mackey 
& Gilmore 2003), in the Fornax Cluster (McLaughlin & van der 
Marel 2005), in irregular galaxies (Larsen et al. 2004) and in inter- 
acting galaxies (Bastian et al. 2006). 



4.3 Star cluster models 

We set up model SCs, ranging in mass from 2 x 10^ M© to 2 x 
lO'' Mq, using an isotropic distribution function (DF): f{x,v) = 
J-{E). The specific form we choose is a lowered polytrope DF 

fix, v) oc [-2E{x, - [-2E^,a.T~'"^ (7) 

with polytrope index n = 2 in all cases. We produce equilibrium 
models through the iterative procedure described in Debattista & 
Sellwood (2000). We set up five such models, C1-C5. All SC mod- 
els have particles of equal mass (1.1 Mq for C1-C2, 5.0 M© for 



4.4 Bare NCD model 

If direct formation of a NCD via gas inflows precedes the full for- 
mation of a NC, how does the accretion of SCs alter the proper- 
ties of a NCD? In order to explore this, we generated a bare NCD 
model by adiabatically growing, over a period of 0.5 Gyr, a disc 
at the centre of the bulge model. The NCD is exponential of the 
type in Eqn. 5, with a scale-length Rd = 9.5 pc and scale-height 
Zd — O.lRd- The disc, truncated at a radius of 5 Rd, consists of 
2 x 10^ particles each with softening e = 0.13 pc . The final to- 
tal mass of the NCD is 1 x 10® M©. We set the kinematics of the 
grown disc to give constant Zd and Toomre-Q = 1.2, as described 
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Table 2. The merger simulations. N(SC) gives the number of stai' clusters used 
and column SC lists which star cluster from those in Table 1 are used in the 
simulation. Column host shows which host galaxy model is used as initial con- 
ditions. 



Run 


N(SC) 


SC 


host 


Comments 


Ml 


8 


CI 


disc 


SCs at mid-plane of main disc 


M2 


8 


CI 


disc 


6 SCs offset from main disc mid-plane 


M3 


8 


C2 


disc 


SCs at mid-plane of main disc 


Al 


10 


C5 


bulge 


multiple accretion of SCs onto NCS 


A2 


20 


C4 


bulge 


accretion of SCs onto a NCD 


A3 


20 


C4 


bulge 


like A2 with 50% retrograde orbits 



in Debattista & Sellwood (2000). For this we calculated the poten- 
tial using a hybrid polar-grid code (Sellwood 2003). 



4.5 Numerical parameters 

All the simulations in this paper were evolved with PKDGRAV 
(Stadel 2001), an efficient, parallel tree-code. We used an open- 
ing angle 9 = 0.7 in all simulations. PKDGRAV is a multi-stepping 
code, with timesteps refined such that 5t — At/2" < ri{e/a)^^'^, 
where e is the softening and a is the acceleration at a particle's cur- 
rent position. We set = 0.1 in all cases. Simulations A1-A3 used 
base timestep At — 10^ years, whereas all other simulations used 
half this value At — 0.5 x 10''' years. 



5 RESULTS OF THE MERGER SIMULATIONS 

We ran three simulations in which eight massive star clusters were 
allowed to merge at the centre of the disc to form a NC. We use 
SC models CI and C2 in these simulations because these SCs are 
sufficiently massive and concentrated to not evaporate too rapidly 
(McLaughlin & Fall 2008) and to have orbit decay times due to dy- 
namical friction less than 3 Gyr from a radius of 1 kpc (Milosavl- 
jevic 2004). The SCs were placed at radii ranging from 14 to 92 pc 
with velocities between 8 and 13 kms~^, which are ~ 60% of 
the local circular velocity. In run Ml the SCs are all initially in the 
mid-plane; in run M2 instead the SCs are vertically offset from the 
mid-plane by up to 67 pc, with tangential velocities similar to Ml 
and vertical velocities up to 1 kms^^. Finally, run M3 is identical 
to run Ml, but uses SC C2 instead of CI. Details of these simula- 
tions are listed in Table 2. 

SCs merge after 13 — 36 Myr, with the longest time needed in 
M2. One SC failed to merge in run Ml and was excluded from all 
analysis. 



5.1 Structural properties 

We measured mass surface density profiles of the remnant NCs 
viewed face-on and obtained _Rcff by fitting Sersic or King profiles. 
The King profile clearly fits the profiles better and we present re- 
sults only of this fit throughout. The King profile is given by (King 
1962) 



T,{R) = k [x^^/^ - C 



-1/2 



(8) 



with normalisation constant k = Eo |^1 — C ^''^j , 
X{R,R,) = 1 + {R/Rcf and C{R,Rc) = 1 + [R/Rcf, 
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Figure 10. Comparison of the simulated and observed NCs in the mean sur- 
face density within R^g versus total NC mass plane. We compare with the 
observed NCs of eaily-type and late-type spiral galaxies as well as Milky 
Way globular clusters. The initial SCs are shown by the open black squares 
while the remnant NCs are shown by black triangles for multiple mergers 
and black circles for multiple accretions. The tracks of the evolving NC in 
the multiple accretion simulations are indicated by dashed lines. 



where Rc is the core radius and Rt is the tidal radius at which the 
projected density drops to zero. Integration yields the cumulative 
form of the King profile for _R ^ _Rt : 



M (R) = TvRtk 



InX 



Cl/2 



1 X 

- + 



c 



(9) 



which is the mass in projection within a cylinder of radius R. _Reff 
can be approximated by 



Rc 



Rc 



1/2 



(10) 



The merger remnants have Red in the range 4.3 — 4.7 pc with 
mass fractions from 87 — 97% of the total merged SC mass. These 
masses are consistent with those of observed NCs (Cote et al. 2006; 
Geha et al. 2002; Boker et al. 2004; Walcher et al. 2006; Seth et al. 
2006). Figure 10 plots the mean surface density within R^s versus 
the total mass of the merger remnants and compares them to NCs 
in early-type galaxies (Cote et al. 2006), late-type spiral galaxies 
(Carollo et al. 1997, 1998, 2002; Boker et al. 2002; Seth et al. 2006) 
and Milky Way GCs (Walcher et al. 2005). As shown by Capuzzo- 
Dolcetta & Miocchi (2008b), we find that SC mergers produce rem- 
nants which have structural properties in good agreement with ob- 
served NCs. 

As was shown already by Bekki et al. (2004), the merger rem- 
nants can be triaxial. We measured the ellipticities viewed face- 
on and edge-on of the simulated NCs using Eqn. 4, obtaining the 
isophote at _Rcff using the task ELLIPSE in IRAF. The NC in Ml 
is significantly non-axisymmetric, with face-on mean ellipticity 
efo — 0.37. We also measured the 3-D shape using the moment- 
of-inertia tensor as described in Debattista et al. (2008). Figure 1 1 
plots the density axes ratios and triaxiality, T — {a^ — h^) / {a^ — (?) 
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Figure 11. The 3-D shape of the NC in runs Ml (soUd lines), M2 (dashed 
Unes) and M3 (dotted lines). The top panel shows the triaxiaUty T, the sec- 
ond row b/a and the third c/a. 



(Franx et al. 1991), of the remnant NCs. In run Ml the NC is triaxial 
within 2 R^s. Models M2 and M3 explore two ways of producing 
more axissnnmetric NCs. In M2 we start the SCs off the mid-plane. 
This makes the NC oblate, with efo — 0.05. In run M3 instead 
we add 2% IMBHs to the SCs which again results in an oblate NC, 
also with £fo — 0.03. 

We find remnant edge-on ellipticities at 3 -Reff £eo in the 
range of 0.36 — 0.56. These values are consistent with the range 
of observed ellipticities 0.39 - 0.89 (Seth et al. 2006). 

Figure 12 shows B4 for the edge-on view. In the triaxial NC 
of Ml, B4 varies with viewing angle but is always negative, i.e. the 
NC is boxy. The NCs in M2 and M3 are also boxy. The merger of 
SCs cannot produce isophotes as discy as observed in the NC of 
NGC 4244. 



5.2 Remnant kinematics 

Bekki et al. (2004) found that his merger remnants were rotat- 
ing, while Capuzzo-Dolcetta & Miocchi (2008a) found that merger 
remnants are kinematically distinct from the main disc/bulge. 

Figure 13 shows the edge-on line-of-sight kinematics of the 
remnant NCs. They are all clearly strongly rotating. However the 
second moment of the velocity, Vrms, shows that the merger rem- 
nants are so dominated by dispersion at the centre that Vrms is cen- 
trally peaked, contrary to what is seen in the NC of NGC 4244 (Fig- 
ure 4). In the bottom row we show the rotation curve Vc (R), the 
line-of-sight velocity V (x), the Une-of-sight velocity dispersion 
<7 (x) and the root-mean-square velocity Vrms (x) of the merger 
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Figure 12. The B4 parameter for the merger simulations. From top to bot- 
tom these are Ml, M2, M3 and the NC of NGC 4244 for comparison. The 
remnant NCs have < 0, whereas the observed NC in NGC 4244 has 
clearly discy isophotes, i.e. positive B4. 
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Figure 14. Final anisotropy (top) and (bottom) in runs Ml, M2 and 
M3. 



remnants. Velocities in M1-M3 peak at larger radii than those ob- 
served in NGC 4244 and M33 (see Figure 1). 

In Figure 14 the anisotropics = 1 — o-'^/aji (top panel) 
and I3z = 1 — a'i /a\ (bottom panel) show that the remnant NCs 
are all radially biased within 4 R^a. M2 is less radially biased than 
M3, which may seem surprising at first, but the radius of the sphere 
of influence of the IMBH in M3 is less than 1 pc, explaining the 
absence of a tangential bias at ReS. The vertical pressure support 
is generally smallest. The initial vertical energy of the SCs in run 
M2 imparts larger vertical random motions, leading to the smallest 
/3z and smallest edge-on ellipticity eeo- Radial anisotropy has been 
noted in the past as a signature of the merging process (e.g. Burkert 
& Naab 2005; Boumaud et al. 2007; Debattista et al. 2008; Thomas 
et al. 2009). The presence of plausible IMBHs does not alter this 
result. 
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Figure 13. Velocity (top row) and Vrms (middle row) fields within 4 R^g for Ml (left), M2 (middle), and M3 (right). The velocity fields show a large scale 
rotation. In all cases, Vrms is centrally peaked. In the top two rows black contours show log-spaced density while the white contours show the kinematic 
contours corresponding to each panel. The bottom row plots the rotation curve Vc (R) (solid lines), the line-of-sight velocities V (x) (dashed lines), the 
line-of-sight velocity dispersions a {x) (dotted lines) and the root-mean-square velocities Vrms (x) (dashed-dotted lines) along the major axis. 



5.3 Accretion onto Super Star Clusters 

The super star cluster (SSC) found by Komei & McCrady (2009) 
in the nuclear region of NGC 253 seems destined to fall into the 
centre of the galaxy and form the basis of a NC. How would the 
accretion of further SCs alter the structure and kinematics of such 
a seed NC? Observed Milky Way globular clusters are found to 
be nearly isotropic (Gebhardt 1994; Gebhardt et al. 1995). How 
much mass needs to be accreted to appreciably alter the isotropic 
distribution? In run Al we study the accretion of SCs onto a NCS 
by introducing a spherical isotropic SC inside the bulge model. We 
form the NCS by letting a massive star cluster fall to the centre. We 
used SC C3 and started it at 127 pc on a circular orbit allowing it 



to settle at the centre over 65 Myrs, before we start the accretion of 
10 SCs. We use model C5 for the accreted SCs, starting them on 
circular orbits at a distance of 32 pc from the centre. In total the 
mass accreted corresponds to ~ 3x the NCS's initial mass. Each 
accretion is allowed to finish before a new SC is inserted. A single 
accretion on average requires ~ 20 Myrs. Table 2 gives further 
details of this simulation. 

After each accreted SC we measure RcS by fitting a King pro- 
file (Eqn. 8) to the mass surface density profile of the NC. The fi- 
nal remnant has R^fi ~ 3.2 pc and structural properties consistent 
with observed NCs, as shown in Figure 10, which tracks the evolv- 
ing NC. The NC becomes triaxial after it has doubled in mass. The 
final efo — 0.17 and eeo — 0.51; the latter is in the observed 
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Figure 15. B4 for edge-on projections of runs Al, A2 and A3 compared 
with tlie NC in NGC 4244 (bottom). For tlie simulations we measured tlie 
isophotes along the edge-on semi-principal axes. 



range (Seth et al. 2006). The final triaxiality T ~ 0.4. Figure 15 
shows that the remnant NC in Al is discy (B4 > 0). 

Multiple accretion of young SCs allows us to also explore the 
effect of a different M / L ratio for a young accreted SC. For the NC 
in NGC 4244 the structural properties are obtained from J-band ob- 
servations (Seth et al. 2006) and the kinematics from A'-band data 
(Seth et al. 2008b). Therefore we obtain both the M/L in the I- 
band and A'-band for the NCD assuming a single stellar population 
with an age about 70 Myrs and a metallicity of [Fe/H] = —0.4 and 
for the NCS assuming two stellar populations, the first with an age 
of 1 Gyr with the same metallicity and the second of 10 Gyrs and 
a metallicity of [Fe/H] = -1.4 in NGC 4244. Using the stellar 
evolution code of Maraston (1998, 2005) this gives M/L ^ 0.2 
for the NCD and M/L ^ 1.6 for the NCS in the /-band and 
M/L « 0.1 for the NCD and M/L ^ 0.8 for the NCS in the 
K-band. Throughout the rest of this paper, we adopt the M/L val- 
ues in the /-band for the analysis of structural properties and M/L 
values in the /('-band for the kinematics, assuming that stars from 
the last accreted SC are young while the rest of the stars are old, to 
obtain a luminosity-weighting. 

Adopting these M / L values, the final edge-on ellipticity be- 
comes eeo — 0.56. 

We produce a luminosity-weighted density map of the final 
NC and fit two components, an elliptical King and an exponen- 
tial disc profile, as in Seth et al. (2006), to measure the struc- 
tural properties of the NCS and NCD component. The NCS has 
RaH ~ 3.2 pc and a flattening e ~ 0.59 while the NCD has 
zo ~ 15.8 pc and scale length Ra ~ 1.6 pc. The NCD consti- 
tutes 3% of the NC luminosity which is a factor of 6 x smaller than 
in NGC 4244. The NCD is ~ lOx thicker than in NGC 4244. 

The initial SSC is isotropic and has no rotation and remains 
unrotating after falling to the centre. Figure 16 shows the evolution 
of (V/a)^; after the first accretion the merger remnant's rotation 
has already increased to (V/a)^ ~ 0.23. By the end of the sim- 
ulation the mass-weighted (V/a)^ ~ 0.34 (~ 0.41 luminosity- 
weighted). The velocity field of the remnant; seen in Figure 17, 
shows that the NC is rotating. However Vrms is still dominated by 
the velocity dispersion at the centre even when luminosity weight- 
ing. The line-of-sight velocity V peaks at larger radii than observed 
in NGC 4244 and M33. 

The settling of the SC to the centre does not change its dis- 




Figure 16. The evolution of (V/a)^ in run Al within R^g (black lines) 
and 3 R^g (red lines). The dashed lines represent mass- weighted and the 
soUd lines luminosity-weighted measurements as described in the text. 



persion within iicff , since the bulge mass within Rett changes only 
about 0.6%. The evolution of and jS^ is shown in Figure 18. 
The first accretion drives /J^ to negative values which slowly in- 
creases as more mass is accreted. The final j3z peaks within RcB 
and declines further out. 



6 ACCRETION ONTO BARE NCDS 

Here we consider the case of a bare NCD, without an initial NCS, 
accreting multiple young SCs. The initial NCD component has a 
mass of 1 X 10® Mq. We use model C4 to represent the infalling 
SCs, placing them at 63pc from the centre of the NC, since the NCD 
is barely affected by the SC until it is well within this radius. In run 
A2 all SCs are on prograde circular orbits, whereas in A3 half of 
the SCs have retrograde orbits. In total the NCD accretes 20 SCs, 
corresponding to 4x its own mass. All SCs orbit in the plane of the 
NCD without any vertical motions. Each accretion event requires 
20 — 30 Myrs to complete; we allow each accretion to finish before 
introducing the next SC. In total these simulations require 1.1 Gyrs. 
Table 2 gives further details of these simulations. 

The remnants both have R^ff — 11.1 pc and structural prop- 
erties consistent with observed NCs as shown in Figure 10, which 
tracks the evolving NC. The surface density increases after each 
accretion event, evolving along the track of observed NCs. When 
we continue to grow the NC in A2 by accreting a further 10 SCs, 
the NC becomes denser than the infalling SC. Figure 19 shows 
the surface density profiles seen edge-on. The top panels show 
mass-weighted maps of the surface density. In the middle row of 
Figure 19 we present luminosity-weighted surface density maps, 
which are more discy than the mass-weighted ones. In the bottom 
panel of Figure 19 we show the surface density profiles. 

Like the merger remnants, the accretion remnants can be tri- 
axial. The NC in run A2 is triaxial with a final efo — 0.22 and 
a T ~ 0.4, whereas run A3 is rounder with efo — 0.05 and 
T ~ 0.1. The final edge-on ellipticity at 3 RcH is eeo — 0.60 
for run A2 and eeo — 0.56 for run A3 and a luminosity-weighted 
ellipticity eeo — 0.73 in A2 and eeo — 0.66 in A3, which is in 
agreement with observed NCs. Accretion leads to the formation of 
NCs with discy isophotes at large radii as can be seen in Figure 15. 
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Figure 17. Kinematic fields for tlie accretion simulations Al (left), A2 (middle) and A3 (right). The top two rows adopt a mass-weighting while the next 
two rows use luminosity-weighting, as described in the text. In the top four rows the black contours show log-spaced density while the white contours show 
the kinematic contours corresponding to each panel. Run A3 shows significant rotation only with luminosity-weighting. The bottom row shows the rotation 
curve Vc (R) (solid lines), the line-of-sight velocities V (x) (dashed fines), the line-of-sight velocity dispersions cr (x) (dotted lines) and the root-mean-square 
velocities Vrms (x) (dashed-dotted lines) along the major axis. ©2011 RAS, MNRAS 000, 1-19 
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R/Reff R/Reff R/Raff 

Figure 18. The anisotropies f}^ (top) and /3z (bottom) after different accretion events for runs Al (left), A2 (middle) and A3 (right). For Al we show the 
initial super star cluster by the black dotted line. For each simulation the red dashed lines show the luminosity weighted anisotropies of the final NC. The lack 
of net angular momentum in ran A3 drives (3^ to negative values, even when luminosity-weighting, whereas in ran A2 f}^ is positive. 
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Figure 19. Projected surface density maps seen edge-on for simulations Al (left), A2 (middle) and A3 (right). The top row shows mass-weighted maps and 
the middle row luminosity-weighted ones. In the bottom panel we show the face-on surface density profiles measured within circular annuli for the NC as stars 
and the corresponding King profile by a dashed line. The bulge is shown by the dotted line and the combined surface density by the solid line. 
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Figure 20. The vertical density profile of simulation A2 for the initial NCD 
(black) and after 20 accreted SCs (red). The solid lines show mass-weighted 
and the dashed lines luminosity-weighted profiles. We plot the density pro- 
file for the total NC, the last accreted SC and the particles of the initial NCD 
in the left, middle and right panel, respectively. 



As in NGC 4244, they show an increase in B4 towards the centre 
and a decline further out. 

The vertical density profile of the NC in run A2 at Rett is 
shown in Figure 20. The NCD is vertically heated by the accretion 
of SCs. The last accreted SC is distributed in a thinner component 
than the initial NCD. We again fit the luminosity-weighted density 
map with an elliptical King and an exponential disc profile as in 
Seth et al. (2006). This gives a NCS with R^h ~ 9.8 pc and a flat- 
tening of g ~ 0.37 and a NCD with zq ^ 1.7 pc and a scale-length 
Rd ~ 3.3 pc. The NCD accounts for 2% of the NC mass, which 
is ~ 3x smaller than in NGC 4244. The scale-height and scale- 
length of the NCD in run A2 is about the same of that observed 
in NGC 4244. Thus, structurally, the NC is comparable to that in 
NGC 4244, provided that the accreted SC is young. 

The initial NCD is strongly rotating. Accretion reduces the 
rotation of the remaining NC. Figure 21 shows the evolution of 
iy/a)^. The first accretion leads to a large decrease in [V/a)^. 
In run A2 the subsequent accretions induce smaller changes, 
quickly asymptoting to iV/a)^ ~ 0.45 (~ 0.52 luminosity- 
weighted), iy/a)^ increases with radius regardless of which 
weighting is used. Run A3 instead drops to {V/g)^ ~ 0.10 al- 
though a luminosity-weighting gives the appearance of more rota- 
tion, {V/a)^ ~ 0.40, comparable to that in run A2. 

Figure 17 plots the first two moments of the line-of-sight kine- 
matics. The remnant NC in A2 is significantly rotating while that in 
A3 shows rotation only when luminosity-weighted. The V shown 
in the bottom panel of Figure 17 peaks at larger radii compared to 
those in NGC 4244 and M33. Unlike in runs M1-M3 and Al, the 
NC remnant in run A2 and A3 has Vi ms that increases with radius, 
like the NC of NGC 4244. However the Vons profiles in Fig. 17 
shows that in run A2 Vrms increases within RcS- In Fig. 22 we 
show that Vrms develops a central peak within _Rcft after the initial 
NCD has doubled its mass. Thus NGC 4244 cannot have accreted 
more than half of its mass as stars. 
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Figure 21. The evolution of the (V/o-)^, within R^g (black lines), 3 -Rcff 
(red lines) and 5 R^g (green fines) in ran A2. The dashed lines repre- 
sent mass-weighted and the solid lines luminosity-weighted measurements 
adopting the M/L of the NCS and NCD in NGC 4244. 
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Figure 22. The Vrms profile after the NC has grown to 1.9 X (black lines), 
2.6 X (red lines), 3.0 X (green lines) and 4.8 X (blue lines) of the initial 
NCD's mass. The solid lines show the Vrms (a;) of the NC seen side-on 
and the dashed lines seen end-on. 

6.1 The vertical anisotropy 

As usual for rapidly rotating discs, the initial NCD is radially bi- 
ased, with /3z initially large but this decreases with accretion, al- 
though it remains positive. Unlike with the kinematic measure- 
ments, no important differences in /3z and fi^ occur if we weight 
by luminosity as shown in Figure 21. Regardless of whether mass 
or luminosity weighting is applied, fiz peaks within 7?cft and de- 
clines beyond. 

The lAM model of the NC in NGC 4244 has a /J^ ~ -0.2. 
We have shown in model M2 that accretion of SCs with verti- 
cal motions decreases fiz- In run A2 the accreted SCs are in the 
plane of the initial NCD. In further tests we let the NC in A2 ac- 
crete the 20th SC on a polar orbit using SC models C4, C5 and C3 
(Af = 2 X 10^ 6 X 10^ 2 X lO^^M©). In Figure 23 we show the 
anisotropy fi^ and fiz of the remaining NC after these accretions. 
Modest accretion off the plane of the disc drives fiz to negative 
values within i?off • The accretion of SC C3 drives /3z < within 
4 i?cff and it also leads to < within RaS- The accretion of SC 
C5 causes *^ within 7?cff snd ■< only within <C 0.5 .Rcff • 
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R/Reff 

Figure 23. The anisotropies (top) and fiz (bottom) for A2 after the ac- 
cretion of 20 SCs and for the thi'ee test runs if we replace the final accreted 
SC with the indicated SC on a polar orbit. 

Thus the observed negative /3z requires that the NC accretes at least 
~ 10% of its mass directly as stars. 



7 JAM MODELS OF SIMULATED NCS 

To compare the rotation of the real NCs and the simulated ones 
in a fully consistent way, we applied the same JAM approach to 
measure the rotation of the simulated NCs. The MGE models were 
fitted to reconstructed images of all models listed in Table 2 at 
an edge-on orientation, for different viewing directions in the disc 
plane. For each projection we fitted the anisotropy (3^ and M/L of 
the simulated NCs, and we then measured the rotation parameter 
K at the best-fitting M/L). The inferred parameters are given 
in Table 3 and compared with the values measured from the parti- 
cles. The comparison shows that the simple JAM models capture 
the global anisotropy of the simulated NCs of the merger simula- 
tions and gives confidence in the values we extracted from the real 
data. In summary for all simulations, as with the observations, we 
measure a degree of rotation k ~ 1 within 5% for models Ml and 
M2, and within 10% for model M3. Given the complex accretion 
process of the NC it is remarkable for the rotation to be so closely 
linked to the NC shape. This result is similar to what was found for 
real galaxies in Cappellari (2008) and suggests that a general pro- 
cess may be responsible for both observations. For the models of 
the merger simulations we recover a weak anisotropy, with models 
for Ml and M3 anisotropic with /3z « 0.25 and ~ 0.17 respec- 
tively, while model M2 is closer to fully isotropic with j3z ~ 0.04 
within 3 RaS- Contrary to models M1-M3, the JAM models do 
not represent acceptable models for the accretion simulations Al- 
A3. In fact the kinematics predicted from the simulated photometry 
under the assumption of constant M/L and an oblate velocity el- 
lipsoid, is qualitatively quite different from the simulated one. The 
main reason for this discrepancy is due to the strong variation in 
the M/L in the simulated NC, which is not included in the JAM 
models. Although it would be easy to construct JAM models using 
the gravitational potential directly measured from the simulations. 



Table 3. Global cylindrical anisotropy and parameters from the JAM models for 
the simulations and observations. We also present the values measured directly 
from the simulations. 







JAM MGE 


Simulation 


Run 


P. 


K 








Ml 


0.25 


0.94 


0.55 


0.15 


0.38 


M2 


0.04 


0.99 


0.37 


0.08 


0.10 


M3 


0.17 


0.90 


0.56 


0.17 


0.26 


Al 


n.a. 


n.a. 


0.56 


0.05 


0.35 


A2 


n.a. 


n.a. 


0.73 


0.10 


0.31 


A3 


n.a. 


n.a. 


0.66 


-1.92 


0.20 


NGC 4244 


-0.2 


1.06 


0.43 






M33 


0.0 


0.84 


0.28 







this cannot be easily done from the observations. However in all ac- 
cretion simulations the global rotation can robustly be determined 
by JAM. 

We perform an independent measure of the degree of rotation 
of the simulated NCs using Eqns. 3 and 4 to place the simulations 
M1-M3 and A1-A3 on the {V/a, e) diagram. For the NCs M1-M3 
the diagram provides results consistent with the JAM models, for 
both simulated and observed NCs. In particular both models MI 
and M3 are in a location of the diagram which indicates signifi- 
cant anisotropy, while model M2 and the observed NC of M33 and 
NGC 4244 are close to the isotropic line (with a typical uncertainty 
of ~ 0.1). The measured values are shown in Figure 5. To first or- 
der the diagram shows that both the simulated clusters and the real 
ones are rapidly rotating. The simulated NCs in runs Ml and M3 
are flatter than M2, but they have a similar (V/a)^. The NCs in the 
accretion simulations Al-3 are similar flat than the merger simula- 
tions MI and M3. They have comparable (V/u)^, although in A3 
only the last accreted SC causes the rotation. All these results are 
consistent with the finding of the JAM models. 

8 DISCUSSION & CONCLUSIONS 
8.1 In situ formation versus accretion 

We have examined in detail the formation of nuclear clusters (NCs) 
via the mergers of star clusters (SCs). This has been proposed as an 
important avenue for NC formation (Tremaine et al. 1975; Mioc- 
chi et al. 2006; Capuzzo-Dolcetta & Miocchi 2008a; Agarwal & 
Milosavljevic 2011). The main support for this mechanism comes 
from the similarity of scaling relations between SCs and NCs (Lotz 
et al. 2001; Walcher et al. 2005). In agreement with previous stud- 
ies (Bekki et al. 2004; Capuzzo-Dolcetta & Miocchi 2008a,b) we 
find that such scaling relations are preserved after mergers . 

As with previous studies (Bekki et al. 2004), the merger of SCs 
was found to produce triaxial NCs, but we showed that axisymme- 
try can result from the presence of intermediate mass black holes 
(IMBH) or sufficient vertical motions. In the only observed galaxy 
where the face-on shape can be determined, M33, we showed that 
the NC is most likely axisymmetric. When the simulated NCs are 
viewed edge-on, mergers produce boxy NCs, unless the merger of 
SCs occurs onto a pre-existing super star cluster or a pre-existing 
nuclear cluster disc at the centre. The flattening is in the range of 
observed NCs in edge-on galaxies (Seth et al. 2006). 

Our simulations indicate that a NC formed via merging of 
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cored SCs leads to a cored NC (Dehnen 2005). During the merger, 
the NC density and mass increases and the NC evolves along 
the track defined by the observed density-mass relation for NCs. 
However the increase in density saturates, when the psc (re) < 
3pNC (re) (Eqn. 8.92, Binney & Tremaine 2008), where r,: is the 
core radius, because the infalling SCs will be disrupted in the outer 
NC, leading to growth in mass and size but not in the mean den- 
sity of the NC, similar to the evolution in mass and size of elliptical 
galaxies due to minor mergers (Ferreras et al. 2009). As seen in Fig- 
ure 10, observed NCs and Milky Way globular clusters (GCs) over- 
lap in the range 10^ < Se < 10* Mopc"^. Some observed NCs 
are denser than the present-day GC population found in the Milky 
Way. Thus if the main formation mechanism of NCs is the accretion 
of SCs, NCs have to form by the merger of denser and more mas- 
sive SCs than those in the Galaxy. Studies of young massive SCs 
in interacting galaxies (Whitmore & Schweizer 1995; McCrady & 
Graham 2007) show that they have similar masses (10® - 10^ Mq) 
and similar sizes as the Milky Way GCs (Bastian et al. 2006, and 
references therein). Due to the smaller infall times of massive SCs 
(Milosavljevic 2004; Neumayer et al. 2011), the mergers of these 
massive SCs is more likely and would explain why NCs are denser 
than the present day GC population. 

Observed NCs contain thin, blue discs of young stars (< 100 
Myrs). In NGC 4244, the mass of this disc is about 5% of the total 
NC mass - if this is typical for the lifetime of a galaxy then over 
a Hubble time dissipation and star formation is sufficient to build 
the NC (Seth et al. 2006). On the other hand we have shown that 
even if the accreted SCs confer no net angular momentum to the 
NC, because the last accreted young SC dominates the luminos- 
ity, the apparent rotation can be quite large. We found that mergers 
can produce rapidly rotating NCs, having (V/u)^ values as large as 
those observed (Seth et al. 2008b). However, the second moment of 
the line-of-sight velocity distribution Vrms, is centrally peaked, un- 
like in the observations. It is only if we introduce a rapidly rotating 
((V/cr)g ~ 2.0) nuclear disc at the centre of the initial system that 
the subsequent evolution produced by infalling SCs is able to qual- 
itatively reproduce the observed Vrms field as well as the observed 
isophotal shape, degree of rotation and mass-density relation, pro- 
vided no more than half of the NCs mass is accreted. Thus a pure 
merger origin of NCs can be excluded. 

Our JAM model of the NC in NGC 4244 has a negative 
f5z = —0.2. Because Cappellari (2008) found that a decrease in 
inclination leads to an increase in fiz, we tested the effect of in- 
clination on Pz for the JAM model of NGC 4244. For inclinations 
> 75°, we found fiz continues to be negative at a 3a level. Thus 
Xhe Pz < in NGC 4244 cannot be a projection effect and must 
be intrinsic. However, a NC disc formed out of gas cooling would 
have I3z > 0. For example, the JAM model of the NC in NGC 404 
has Pz ~ 0.5 (Seth et al. 2010); in this case the observed burst 
of star formation ~ 1 Gyr ago could easily have fed its NC and 
given rise to its positive /Sz. By scattering box orbits, central black 
holes can lower I3z (e.g. Merritt & Quinlan 1998). Whereas Seth 
et al. (2008a) find indirect evidence of supermassive black holes 
(SMBHs), with masses ranging from 10-100% in - 10% of NCs, 
in NGC 4244 the radius of influence for the largest possible inter- 
mediate mass black hole is < 1 pc, while the NC has Res = 5 pc. 
Therefore the presence of a intermediate mass black holes in the 
centre of NGC 4244 cannot explain its !3z < 0. The only way we 
were able to produce /3z < was through the accretion of SCs on 
highly inclined orbits. Thus, at least ~ 10% of the mass of the NC 
of NGC 4244 needs to be accreted as SCs in order to obtain 0z < 0, 



which constitutes a lower limit on the amount of mass accreted in 
the form of star clusters. 



8.2 Nuclear Cluster Formation in Dwarf Elliptical Galaxies 

Our simulations can also make testable predictions relating to the 
formation of nuclei in dE galaxies. Lotz et al. (2001) suggested that 
NCs in dE galaxies could have formed by the merger of GCs. They 
foimd a depletion of the most massive GCs relative to the less mas- 
sive ones in the inner region of galaxies. They interpreted this as 
being due to shorter dynamical friction infall timescales for mas- 
sive GCs than for the lower mass ones. In the ACS Virgo Cluster 
Survey, Cote et al. (2006) fotmd that, on average, NCs in dE galax- 
ies are 3.5 mag brighter than the mean GC. So if NCs form via the 
merger of GCs, on average about 25 GCs need to merge. Similar 
to what we found for NCs in late-type spiral galaxies (Figure 10), 
Cote et al. (2006, see their Figure 18) found that the NCs are denser 
than the mean density of the GC population in dE galaxies. They 
also found that the fraction of red GCs and the (g — z)'^j^ colour 
of NCs increases with the host galaxy luminosity (Peng et al. 2006; 
Cote et al. 2006). Using Monte Carlo simulations, and assuming 
that these NCs formed via mergers, Cote et al. (2006) found that the 
resulting scaling relation of the NCs colours is less steep than those 
observed. From spectroscopic data Paudel et al. (2010) found that 
the metallicity of NCs in a sample of dE galaxies correlates with the 
luminosity of the host galaxy. They also found that the median dif- 
ference in age between the NC and the galactic main body is about 
~ 3.5 Gyrs and that the difference is more prominent in dEs with 
discy isophotal shapes. This implies that the formation of NCs in 
dEs might be enhanced by the accretion of gas. In contrast, Paudel 
et al. (2010) found fairly old and metal-poor NCs in very faint dEs, 
resembling the properties of the GC population. This suggests that 
NCs in faint dEs might have formed by different processes than the 
NCs in brighter dEs. 

In summary both the accretion of gas with in situ star for- 
mation and the merger of GCs could be at work to form NCs in dE 
galaxies. If NCs in dEs form via merger of GCs, our simulations in- 
dicate that they will have boxy shapes, have centrally peaked Vrms 
and be radially biased. On the other hand, if the main formation 
mechanism is dissipation, NCs will have discy isophotes and no 
centrally peaked Vrms- At present, no available observational data 
exists which is able to test these predictions. 

8.3 The McMO-o"e relation 

The star formation histories in NCs of late-type spiral galaxies are 
extended, with the yovmgest population of stars less than 100 Myrs 
old (Walcher et al. 2006; Rossa et al. 2006; Seth et al. 2010). NCs 
appear to be offset from the M, — Ue relation of SMBHs (Ferrarese 
et al. 2006a; Wehner & Harris 2006). NCs have the same slope, 
but, for a given velocity dispersion a^, are lOx more massive than 
SMBHs. McLaughlin et al. (2006) found that this offset can be ex- 
plained by feedback from stars and supernovae. Our simulations 
indicate that at least 10% of the mass of NCs needs to be accreted 
by SCs and that at least 50% of its mass needs to be accreted as 
gas. Therefore a fraction of the star formation could occur outside 
the NC. However, the accreted SCs have to be young and therefore 
still formed within the central region. SCs with masses in a range 
of lO'^ - 10'^ M© have to be formed within 60 - 200 pc, otherwise 
their infall times are longer than 100 Myrs. Therefore, even if the 
NC has accreted half of its mass in yotmg SCs, stellar feedback oc- 
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curs within tlie central region of tlie galaxy and therefore remains a 
plausible explanation for the Mnc — fe- 

8.4 Summary 

We have studied the formation and evolution of stellar nuclear clus- 
ters using A'^-body simulations. Our main conclusions can be sum- 
marised as follows: 

• We find no evidence of non-axisymmetry in the nuclear cluster 
of M33. Its PA is consistent with that of its main disc and its appar- 
ent ellipticity is consistent with a vertical flattening q = 0.7, which 
is the average observed in the NCs of edge-on late-type galaxies 
(Seth et al. 2006). There is also only a small misalignment between 
the photometric and kinematic PAs. While this is the only galaxy 
where this measurement can be done at present, it suggests that 
NCs are generally axisymmetric. 

• The NC in NGC 4244 is nearly isotropic, with (3^ = -0.2 ± 
0.1, and is rapidly rotating. It has a mass of (1.1±0.2) X lO'' Mq; 
if an IMBH is present its mass is less than 1% that of the NC. It is 
not possible from these models to distinguish whether the rotation 
is present throughout the NC or is restricted to the NCD. 

• The mergers of SCs produce NCs with density and sizes con- 
sistent with observations, evolving along the track of observed 
NCs. Remnant NCs can be axisymmetric. Multiple accretions of 
yoimg SCs onto a pre-existing nuclear cluster spheroid can also 
produce discy isophotal shapes. Mergers can lead to rapidly rotat- 
ing NCs, as observed. They have 0z ^ 0.4, where initial vertical 
motions induce the smallest value. However these NCs have cen- 
trally peaked Vrms unlike observed NCs in late-type discs. 

• Accretion of young SCs onto a bare NCD produces NCs with 
densities, sizes and ellipticities comparable to those observed. They 
show discy isophotes and have rotations comparable to those in the 
NCs of NGC 4244 and M33. Vrms is dominated by dispersion if 
the accreted mass is greater than the initial mass of the NCD. The 
formation of the NC in NGC 4244 therefore requires at least 50% 
of the mass to be accreted as gas to match the observations. 

• The negative in the NC of NGC 4244 requires at least ~ 
10% it's total mass to have been accreted from SCs. 

• We caution that even if the accreted SCs do not impart any net 
angular momentum, the last accreted SC can dominate the appar- 
ent rotation when luminosity-weighting and could be similar to the 
C^/'^)e Vrms fields of obscrved NCs. 

• With the results presented in this paper, the simulations are 
now ahead of the observations with predictions of detailed observ- 
ables that can be used to constrain the formation scenarios better. 
Integral-field observations of the kinematics of more NC are essen- 
tial for further progress and we hope these will be obtained in the 
near future. 
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